Resummed Perturbation Theory of Galaxy Clustering 
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The relationship between observed tracers such as galaxies and the underlying dark matter dis- 
tribution is crucial in extracting cosmological information. As the linear bias model breaks down at 
quasi-linear scales, the standard perturbative approach of the nonlinear Eulerian bias model (EBM) 
is not accurate enough in describing galaxy clustering. In this paper, we discuss such a model in 
the context of resummed perturbation theory, and further generalize it to incorporate the subse- 
quent gravitational evolution by combining with a Lagrangian description of galaxies' motion. The 
multipoint propagators we constructed for such model also exhibit exponential damping similar to 
their dark matter counterparts, therefore the convergence property of statistics built upon these 
quantities is improved. This is achieved by applying both Eulerian and Lagrangian resummation 
techniques of dark matter field developed in recent years. As inherited from the Lagrangian descrip- 
tion of galaxy density evolution, our approach automatically incorporates the non-locality induced 
by gravitational evolution after the formation of the tracer, and also allows us to include a continu- 
ous galaxy formation history by temporally weighted-averaging relevant quantities with the galaxy 
formation rate. 



I. INTRODUCTION 



Large-scale structure surveys provide a wealth of cos- 
mological information in probing dark energy, modified 
gravity, neutrino masses and the physics of early uni- 
verse. As the statistical uncertainty decreases dramati- 
cally in next generation surveys, it requires us to achieve 
an unprecedented level of accuracy in theoretically pre- 
dicting the statistics of the observed clustering pattern. 
Since nonlinear dynamics plays a important role in un- 
derstanding the observed data in relevant regime, it has 
attracted lots of interest in developing the nonlinear per- 
turbation theory of the dark matter beyond the stan- 
dard approach. Many sophisticated methods have been 
proposed, including the renormalized perturbation the- 
ory (RPT) [IH3; the Lagrangian resummation theory 
0], the closure theory 0], and the time renormalization 
group theory (TRG) 0. However, the relationship be- 
tween various observable tracers (galaxies, quasars, Lya 
forest, clusters and HI galaxies etc.) and the underlying 
dark matter field may be more challenging. The linear 
bias model breaks down even for present surveys @, H| 
at relatively large scale, the scale-dependence observed 
in the data varies among different galaxy samples as one 
would expect Q, and the empirical fitting formula @ 
widely used in the data analysis is not able to describe 
such differences among samples and eventually leads to 
inconsistent cosmological constraints Q- 

Two different approaches exist in modeling the scale- 
dependent bias of galaxies these days. One is the halo- 
based model which starts from the large scale clustering 
of dark matter halos and then populates with different 
types of galaxies according to the halo occupation dis- 
tribution (HOD) function. The other approach is the 
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perturbative bias model, which packs the complicated 
baryonic physics into a general unknown functional de- 
pending on the dark matter distribution, and character- 
ize it with the first few bias parameters. With a careful 
parameterized and calibrated HOD function, the first ap- 
proach can describe the two-point statistics qualitatively 
well at small scale. However, a viable HOD parameteri- 
zation usually contains lots of free parameters, and heav- 
ily depends on calibration from simulation. Moreover, 
in the quasi-linear regime where most of large-scale sur- 
veys probe, the accuracy of such approach is usually not 
satisfactory. 

On the other hand, the standard perturbative bias 
model is also problematic. In the Eulerian bias model 
(EBM), the observed galaxy number density depends on 
the Eulerian dark matter distribution at the time of ob- 
servation. Therefore, this model breaks down as the dark 
matter perturbations evolve into the nonlinear regime at 
small scales and later times. Since both the large-scale 
amplitude and small-scale shape of the power spectrum 
depend on higher-order bias parameters, such calculation 
would fail within the whole range of scales. Furthermore, 
although adopting some artificial nonlocal bias parame- 
ters won't bring too much intrinsic difficulties in the cal- 
culation, most authors simply assume that the functional 
form of the galaxy density only depends on the matter 
distribution locally, which is not consistent with the non- 
local gravitational evolution of the tracer [Io|, HH ■ 

Since the difficulties of the perturbative convergence 
encountered in the nonlinear bias model is very similar 
to the situation in dark matter field, one would expect 
to generalize any of previous nonstandard matter pertur- 
bation theory to describe the clustering of biased tracer. 
Several works have already been done along this direc- 
tion. In [ljj], a Lagrangian bias model (LBM) which in- 
tended to incorporate the halo bias and redshift distor- 
tion was proposed based on a similar resummation the- 
ory for matter distribution 4]. Then [l3[ also included 
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the halo bias in the time renormalization group methods. 
[lil [l5| approached in a different way by redefining the 
bias parameters. As we will show in this paper, with the 
help of recent development in the renormalized pertur- 
bation theory [J-Q, it is also possible to construct a re- 
summed perturbation theory for nonlinear Eulerian bias 
model, and the perturbative expansion of statistics can 
be rewritten through the so-called T— expansion, first in- 
troduced by [|[ , with the building blocks known as multi- 
point propagators. One advantage of the T— expansion is 
that for Gaussian initial conditions, all contributions are 
positive and centralized within a limited range of scale. 
As we will see later, such an expansion is quite general 
and applicable to many circumstances as long as the mul- 
tipoint propagators can be well estimated. Therefore in 
the following, this paper will mainly concentrate on con- 
structing the multipoint propagators for various objects. 
The validity of such technique has already been verified 
in calculating the power spectrum of logarithmic map- 
ping of the density field [lg|. It have been shown that 
even for the slow converging function such as logarithm, 
the calculation agrees with the simulation quite well. 

We further generalize the Eulerian bias model by sepa- 
rating the complicated nonlinear physical processes dur- 
ing the galaxy formation and the subsequent gravita- 
tional evolution until been observed by a survey. As we 
will show, this can be achieved by combining the Eu- 
lerian bias model and Lagrangian description of galaxy 
motion introduced in [12|]. Like EBM, we assume the 
galaxy distribution created at some time is characterized 
by a general unknown functional, and the galaxy number 
density observed at a Eulerian position later can be nat- 
urally described by the Lagrangian version of continuity 
equation. Consequently, this approach also enables us to 
incorporate continuous galaxy formation by temporally 
averaging the distribution of newborn galaxies weighted 
with the galaxy formation rate. In such a framework, 
while the complex galaxy formation physics is simplified 
by a few bias parameters, we are able to capture the fol- 
lowing gravitational evolution after its creation, which 
would in principal induce non-locality and can not be ac- 
curately described by a local Eulerian bias model. As in- 
herited from the Lagrangian bias model and our resumed 
calculation of Eulerian bias model, this model can also 
be expressed in a resummed manner. 

In section II, we briefly review both Eulerian and La- 
grangian perturbation theory of dark matter field, includ- 
ing both standard approach and resummed theory. We 
will also introduce the concept of the multipoint propa- 
gator in the context of Lagrangian perturbation theory, 
and give a explicit calculation in Appendix A. In section 
III, we first review the standard Eulerian and Lagrangian 
bias model, and then introduce our continuous galaxy 
formation model. We construct the resummed Eulerian 
bias model in Section IV, and the generalize it in Section 
V. Finally, we conclude in Section VI. 



II. PERTURBATION THEORY OF DARK 
MATTER 

In this section, we will briefly review both Eulerian 
perturbation theory (EPT) and the Lagrangian pertur- 
bation theory (LPT) of the dark matter field, including 
both the standard approach and resummed calculation. 



A. Dynamical Equation of Motion 

The gravitational dynamics of a pressureless fluid be- 
fore shell crossing is governed by the continuity, Euler, 
and Poisson equations. 



&S(x,t) 

dr 
9u(x, r) 

dr 



+ V- [ (l+<5(x,r)) u(x,r) ] =0, 

+ H(t)u(x, t) + u(x, t) ■ Vu(x, r) 

= -V$at(x,t) 



V 2 $ jv (x,t) = -n m (r)n 2 (r)S(^r). 



(1) 



Here, 7-L = d\na(r)/dhiT is the Hubble expansion rate, 
a(r) is the scale factor, f2 m (r) is the ratio of matter den- 
sity to critical density, and <fr at ( x , t) is the Newtonian 
potential. 

Following [H,0|, the equation of motion in Fourier space 
can be written in a compact form by defining the two- 
component variable 



H a (k,77) = (5(k,77), -6(k lV )/H), 



(2) 



where 9 denotes the divergence of peculiar velocity V • u. 
The index a £ {1,2} stands for matter density or ve- 
locity variables, respectively. In the following, we will 
interchangeably use Si, S m and S m for matter overden- 
sity, while using both S g and 5 g for galaxy perturbation. 
The equation of motion then reads 

d v S a (k, 77) + f2 af ,H 6 (k,?7) 

= 7abc(k, ki,k 2 )Hb(ki,7])H c (k 2 ,?y), (3) 

with the convention that repeated Fourier arguments are 
integrated over. The time 77 = In a(r) in a Einstein-de 
Sitter (EdS) universe. Here the constant matrix 



-1 
-3/2 1/2 



(4) 



derived for an EdS universe, is still applicable in other 
cosmologies with negligible corrections to the coefficients, 
using r/ = In D(t) (with D the linear growth factor). 
The symmetrized vertex matrix Jabc is given by 



7222(k,ki,k 2 ) = 5 D {k-ki - k 2 



7i 2 i(k,ki,k 2 ) = 6o(k-ki-k 2 ) 



|k! + k 2 | 2 (k! ■ k 2 

^2 



(ki + k 2 ) ■ ki 

2kf 



(5) 
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7m(k,k 1) k 2 ) = 7i 2 i(k,k 2 ,ki), and 7 = otherwise. 

Then the formal integral solution to Eq. ([3]) can be 
derived as 



3a(k,77) 



9ab{ri)<l>b(k) + / drf g a b(v - v') 
Jo 

xj bcd (k, ki, k 2 )S c (k 1 , 57')a d (k 2 , rjf) (6) 



where </> a (k) denotes the initial condition </> a (k) = 
S a (k, rj = 0), and the linear propagator g a b(v) is given 
by 



Jab(v) 





'3 2" 




-2 2" 


y 


3 2 


5 


3 -3 



(7) 



In the following, we adopt growing-mode initial condi- 
tions a (k) = <5 (k)u a , with u a = (1, 1). 

In Lagrangian perspective, the dynamic of a fluid ele- 
ment at initial Lagrangian position q is entirely described 
by the final Eulerian position x, or equivalently the dis- 
placement field \I/(q) = x q, which is governed by the 
equation of motion, 



(8) 



where V is still the derivative with respect to Eulerian 
coordinates x, and the gravitational potential $n is de- 
termined by Poisson equation as shown in Eq. ([1}. The 
density contrast is then related to ^(q) via mass conser- 
vation 



i + Mx) = J- 1 



det [6ij +*ij(q)] 



(9) 



where J is the Jacobian of the transformation between 
Eulerian and Lagrangian space, and ^>i j = d^i/dqj. To 
linear order, Eq. © is then solved by the Zel'dovich ap- 
proximation Vg-^^(q) = — D(t)So(<1), with Lagrangian 
derivative V Q . 



B. Eulerian Perturbation Theory 

A perturbative solution to Eq. © could be obtained 
by expanding in terms of initial fields 

00 

S (M) = X>£°(M) 

n=l 

BW(k,7?) = |rfV..nfe(k- qi ... n )J-W.. b „(qi,---, 

• • • , q«; v)<t>h (qi) • • • 4>b n (q«) (10) 

where d 3 qi... n is short for d 3 qi • • • d 3 q„, and qi...„ de- 
notes qi + • • • + q n . The kernels F^ are fully symmetric 
functions of the wave vectors. As shown in [2j, they can 



be obtained in terms of g ab and 7 aDC recursively. 
J"a n) (ki,-- - ,k n ;r?)5 (k-kx... n ) = 
ELi Jo drfgabin ~ i) T6cd(k, ki... m , k m+1 ... n ) 



x J~c (k]..., 



(11) 

symmetrized 



For n = 1, F a {rj) — g a b(r]) u b- It should be noted that in 
this formalism, the kernel depends on the time 77, since it 
includes subleading terms in e v [l| . If one only considers 
the fastest-growing mode, F a n ^ {rj) equals the well-known 
PT kernel D n (rj){F^ ,G^}, where D(rj) is the linear 
growth factor. 



-®- 




FIG. 1: Diagrammatic elements of Eulerian perturbation the- 
ory, left: perturbative expansion (Eq. [10|; right: initial power 
spectrum Po(fc). 

The power spectrum P m (k,r)) of the matter density 
perturbation is defined as 

P m (M)Mk + k') - (SmCM) Sm(k',»7)> (12) 

For Gaussian initial conditions, all the statistical infor- 
mation is encoded in the initial power spectrum 



(0 Q (k)<Mk')> = ^(k + k , )P ofe (fc) ) 



(13) 



where P ab (k) = u a u b P {k), with P (k) = (<5 (k)o\,(k')). 
Diagrammatically, as shown in the right-hand diagram 
of Fig. ([U , the ensemble average is obtained by gluing 
two open circles together to form the symbol ®, which 
represents the initial power spectrum P ab {k). 

As the density fluctuations evolve into the non-linear 
regime at later times, the validity of standard perturba- 
tion theory breaks down, loop contributions become ill- 
behaved and the convergence of perturbation series gets 
out of control. As one of several different approaches be- 
yond SPT been proposed recently, we will briefly review 
the renormalized perturbation theory (RPT) introduced 
by [H-Q, and then the Lagrangian resummation theory 
1 1| in the next subsection. 

The crucial step of renormalized perturbation theory 
is to define the generalized growth factor, known as the 
propagator ab {k) 



r£U(M)Mk-k') = 



_ / <5E (k,7?) 



(14) 



o— « — = — « — ■ — * — + o y < — + <> < < — + 



FIG. 2: The nonlinear propagator F„ (fc, if) has an infinite number of loop contributions. 



which effectively describes the time evolution of individ- 
ual Fourier modes when non-linear mode-coupling is in- 
cluded. Here 8 denotes the functional derivative. The 
propagator measures the dependence of a non-linearly 
evolved Fourier mode H (k, rj) on its initial state </>t>(k) 
on average. Intuitively, one should expect G a b to decay 
to zero at small scales since non-linear mode-coupling has 
erased all the information from the initial state at that 
k. Indeed.with the help of the Feynman diagrams intro- 
duced in |l|, the dominant contribution can be summed 
up explicitly in the large-k limit, and giving the exponen- 
tial decay G a b(k) « exp(— k 2 D 2 a%/2). In this paper, the 
non-linear propagator is diagrammatically presented as a 
grey circle with an incoming branch. It is a summation 
of infinite number of loop contributions, as illustrated in 
Fig. ©. 

Given above, [lj was able to rewrite the nonlinear mat- 
ter density power spectrum as a sum of two contributions 



known as the T— expansion 



P m (k,v) = G 2 (k, v )P (k) +P M c(k,r,), 



(15) 



where G is the density propagator G(k 7 r]) = GuUb, and 
PMc{k 7 i]) is the mode-coupling term. Therefore, the 
non-linear power spectrum at any k is composed of two 
parts. One is proportional to the initial power spectrum 
at the same k; the other comes from mode-coupling of 
other k'. As G(k) decays at small scales, more and more 
power comes from the mode-coupling contribution. [3] 
showed that these complicated mode-coupling contribu- 
tions can be expressed as a summation of multi-point 
propagators, defined as 



1 



,k„;^)(5_D(k-ki... n ) 
8 n 3 a (k, V ) 



bl (ki) 



,(k„) 



(16) 



which is nothing but a generalization of two-point prop- 
agator G a b- Similarly, the dominant part of multi-point 
propagators can also be summed and decay into the non- 
linear regime at the same rate as two-point propagator. 
Q showed that a simple approximation which general- 
izes the fc-dependence of two-point propagators agrees 
with the data with acceptable accuracy, 



rW(k Xj -.. ,k n ; V ) = r i ^ trec )(k 1 ,..- ,k n ;7?) 

^rn ( |k]_... n |) 



t(1, tree) 



(|ki 



(17) 



For gaussian initial condition, the nonlinear power spec- 
trum can then be expressed with the help of Ea.dTI 



P m ab(k,rj) = y^n! J d 3 qi... n 8 D (k - qi...„) 



r ma(qi> - ,q»;f?) r^ 6 (qi,-- - ,q„;r?) 



P (<Zl)---P (?n). 



(18) 



Note that Eq. (|18[) describes nothing but an alter- 
native way of taking ensemble averages, or diagram- 
matic speaking, gluing initial states. Instead of gluing 
two density fields order by order, one can first construct 
the objects by gluing initial states of individual density 
fields with n incoming branches, known as an (n + 1)- 
point propagator, and the final non-linear power spec- 
trum is then obtained by gluing two propagators to- 
gether. Therefore, this T— expansion is quite general and 
applicable to the ensemble average of any two arbitrary 
statistical fields (£1X2), which are functionals of some 
other field x$ 

Xi [x (lL n ),X\ = J2 I d ^i-n ^ (n) (ki,--- ,k„;A) 

x [x (ki)---x {K)], (19) 

as long as the building blocks are well estimated. Here 
A symbolize some parameters, and X( n > are symmetric 
on their arguments. In the case that xq is Gaussian, 
the Wick's theorem ensures that the ensemble average 
of an odd number of fields vanishes (xq(1) ■ ■ ■ xo(2n + 
1)) = and the average of an even number of fields 
(xq(1) ■ ■ ■ xo(2n)) can be decomposed as summation of 
two-point correlations 

(x (l)---x (2n)>= ]T H 

(x (i)x (j)). (20) 

all pairs ij 

Following the idea of Q, there are three types of two- 
point correlations (xo(i)xo(j)) involved in calculating 
(xi(X)x2(X'))'- pairs within x\ or X2 themselves and pairs 
connecting them. Therefore, one is always able to re- 
group Eq. (|2"0)) in the following way 



(1)... 10(270)= £ J2 EiMDMm)) 



in between 



x\ internal Irn 



e n^>o(™')> n <*o(osocj)> 

ij in between 



internal V 



(21) 
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Combining above equations (TT9ll20l and[2T j) together, one 
is able to reformulate (xi(A)x2(A')) as a sum of two-point 
cross correlations of Xo belonging to x\ and X2 respec- 
tively with some pre-summed nonlinear quantity involv- 
ing internal pairs £ intornal Yli m ( x o(l)xo(m)) , which can 
be formally defined as the average of functional derivative 



n\ \ (5xo(l) • • • 5xo(n) 



(22) 



Therefore, as long as r\ is well estimated, one is always 
able to carry out the T— expansion in such situation. In 
the following, we will demonstrate that with the help of 
the density propagator of cold dark matter (Eq. ITB"]) . 
many similar quantities of this kind can be calculated. 
And we will extensively utilize such expansion in the fol- 
lowing of the paper. 




For the matter distribution, one usually assumes that the 
initial field is sufficiently uniform 5q ~ 0, so that above 
equation (Eq. |24|) can be further simplified. In Fourier 
space, 



6 m (k) = J d 3 q e- lk i 

= H)" 



= -ik-*(q) 



1 



<i 3 pi...„ [k- *(pi) 
x S D (k-p 1 ... n ] 



[k"*(Pn)] 

(25) 



Similarly, the above equation can also be represented by 
diagrams as shown in Fig. ([3]). Here, wavy lines denote 
the displacement field and the thin solid line together 
with grey open circle stands for the linear density per- 
turbation, i.e. Eq. (|25|) . 



o 





FIG. 3: Diagrammatic representation of nonlinear matter 
density perturbation in Lagrangian perturbation theory, sim- 
ilar to the one introduced in [J]. The wavy line denotes the 
displacement field ^f, and the thin solid lines together with 
grey open circles stand for the linear density perturbations 
S . 



C. Lagrangian Perturbation Theory 

As discussed in the beginning of this section, the dy- 
namics of the system in Lagrangian picture is fully char- 
acterized by the displacement field ^(q, rf). Combining 
the equation of motion (Eq. [5J and the Poisson equation, 
the displacement field can be solved perturbatively. In 
Fourier space, the n— th order perturbation *!>(") can be 
generally expressed as 



*W(k) = —r- / d 3 Pl ... n I»( Pl ... n )<Jo(pi 
n\ J 



•<*o(Pn), 



(23) 



where LA") is the perturbation kernel. 

Due to the matter conservation, the Eulerian density 
contrast 5 m (x) can be expressed as 

1 + S m (x) = y"d 3 q[l + (5o(q)](5 I? (x-q-*(q)). (24) 



FIG. 4: Two-point propagator of matter density r„'(k) in 
the context of Lagrangian perturbation theory up to one-loop 
order. The grey ellipse shown in each contribution symbol- 
izes the quantity of Eq (|26|l . Note here we have omitted the 
summation sign over the displacement field VE* presented in 

Starting from Eq. (|25[) . Q performed a resummation of 
the power spectrum of the density contrast, known as La- 
grangian Resummed Perturbation Theory. Although not 
been explicitly discussed in Q , it is also possible to define 
the same propagator within this framework. First of all, 
as pointed out in [l(| for biased tracers 1 , the following 
quantity can be resummed 

°2, (—i) n f 

/ d3 Pl-« fci •••*<»] (*ii(Pl) (Pn)) 



E 



d 3 P k- *( P ) 



ik-*(0) 



>, (26) 



where \P(0) = J d 3 p \P(p) is the displacement field at the 
origin. It is represented as a grey ellipse in Fig. (|4]), and 



1 Setting bi = 1 and b n = otherwise in Eq. (71) of |lCj . one 
recovers the following equation. 
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can be further simplified as 0, [T(| 
-ik.*(o)\ = n (°)(k) =exp 



-ik-*(Q) 



exp 



exp 



~ (-l)"(|^(Q)p") n 
^ (2n + l)(2n)! 

n=l v f \ / 

k 2 D 2 



(27) 



with the variance of displacement field along one direc- 
tion D 2 a% = (|*(0)| 2 )/3 = Jd 3 qD 2 P (q)/(3q 2 ). Here, 
we have used the cumulant theorem 



[5Dn>0 3nX n ] \ _ 



exp 



;[£„>0 j" 3 ^] _ 1 



(28) 



where (• • • ) c denotes the connected part of the ensemble 
average. As shown in (Toj . it leads us to a partial resum- 
mation of the propagator as represented in Fig. Q up to 
the one-loop order. Again, this is just a simplified version 
of Fig. (17) in [10] for the Lagrangian bias model. The 
grey ellipse represents all possible graphs attached with 
various number of wavy lines, i.e. Eq. (1261) . Compared 
to the diagram introduced in [Toj . we have omitted the 
summation sign over external wavy lines. Mathemati- 
cally, such diagrammatical series is equivalent to 



r£>(k) 



-ik-*(0) 



d 3 p 



-ik- 



7T*(P) 



IT:ni,k ^(^M expHk '* (0)] ) c 



(29) 



Note that we have used an alternative version of the cu- 
mulant theorem in Eq. ([28)) . 

(e jlXl x 2 ) = exp [(e^ 1 - l)J {{e 3lXl x 2 ) c } . (30) 



which is obtained by taking the derivative of Eq. (|28|) 
respect to ji and then setting it to zero. We further cal- 
culate these contributions explicitly in the Appendix A. 
The result is not the same as the form originally proposed 
by ■ However, the matching prescription between high- 
k and one-loop calculation in Q seems more or less ad 
hoc. Therefore, an alternative way of constructing the 
propagators was proposed by [2(J, to which our calcu- 
lation do agree. See [2(| for the discussion of different 
versions of calculating multipoint propagators. 



III. THE NON-LINEAR BIAS MODEL 



A. Eulerian Bias Model 

In Eulerian bias model, the galaxy number density per- 
turbation at time r\ is considered as some general func- 
tional F[<5 m ] of the underlying matter density distribu- 
tion characterized by its Taylor coefficient b n 

Mx, v) = f [Mx, n)] = y\— y Pm(x, r,)T , (31) 




Pl,l 

pi ,2 

/'!.:( 



P«,l 



FIG. 5: Diagrammatic representation of fe(k, rf) in Eulerian 
bias model. 



where the nonlinear mass density field 5 m (x) is solved 
by perturbation theory mentioned in the last section. As 
discussed by [13], this functional is generally non-local. 
Therefore in Fourier space, bias coefficient b n depends on 
wavenumbers 

°° 1 f 

fe(k,r?) = ^2 — / d 3 ki- n S D (k - ki...„)6„(ki, • • • ,k„) 



x [S m (ki ,77) • • -5 ri)\- 



(32) 



In the local bias model, b n then becomes constant 
times a smoothing window function W(k,R) in Fourier 
space. In this model, the linear bias relation at large 
scale P(k) « b 2 P m (k) is not only attributed to the first 
term in Eq. (|32|) . but also from the higher orders b n . As 
first pointed out by [13], the effective large scale bias is 
determined by infinite number of higher-point local vari- 
ances of the field. Therefore, similar to the situation in 
the standard perturbation theory of matter density field, 
at later time when mass density variance grows larger, 
a naive perturbative expansion fails. In this case, such 
a breakdown affects both the large scale amplitude and 
the high-k regime of the power spectrum. 

Given the definition of fe(k) in Eq. ([32|) and the dia- 
grammatic representation of the density perturbation in 
Fig. (JTJ) , it is also possible to draw diagrams representing 
fe(k). As depicted in Fig.©, each diagram contains two 
levels of nonlinear interactions. The thick solid n-branch 
tree represents the nonlinear convolution, and their inter- 
action vertex, shown as a solid square, carries the bias co- 
efficient b n /n\, where n is the number of branches. Each 
branch represents one nonlinearly evolved density con- 
trast 5 m (]s. n ,ri) at time 77, which is followed by another 
m thin branches, to express its gravitational nonlinearity 
in terms of initial fields and perturbative kernels J- i - m \ 
At the end of each branch, small open circles represent 
the initial conditions <5o- When expanding 5g(k, 77) to 
different orders, one should also include the information 
about the number of topologically equivalent diagrams, 
which essentially comes from the multinomial coefficient 
of (6 m ) n . The final field fe(k, rf) is then a summation of 
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all possible diagrams. 



B. Lagrangian Bias Model 

To incorporate the halo bias with nonlinear pertur- 
bation theory, fl2j proposed an alternative bias model 
involving the initial Lagrangian density fluctuation S (q) 
instead of final Eulerian density 6 m (x, 77) . In this La- 
grangian bias model, the non-linear objects, typically 
halos, were formed at some initial time with the num- 
ber density is a functional _F[Jo(q)] of underlying linear 
density distribution 8q. Then the Eulerian observable 
<5_l(x, 77) at later time is described by the continuity rela- 
tion 



1 + 5l(x) 



d 3 q [1 + S L (q)} fo[x- q- ¥ L (q)] 



J d 3 qF[6 (q)} fe[x-q-*i(q)], 



(33) 

where ^i(q) is the displacement field of such tracer. In 
Fourier space, 

S L (k,r]) = ( d 3 q e- jk - q \F[S Q (q)} e^*^ - ll .(34) 



In the limits that nonlinear objects trace exactly the un- 
derlying dark matter, i.e. ^ l — ^m, one can further 
proceed to calculate the power spectrum P(k) using La- 
grangian perturbation theory of the matter field. In Fig. 

we reproduced the diagram of this model introduced 
by flp| on the left, where the wavy line represents the 
displacement field *S?l and the n branches of thin solid 
lines with open circle symbolize the nth order of -F[$o]- 

The great advantage of this approach is that it natu- 
rally takes account of the non-locality induced by gravita- 
tional evolution. As shown by [Toj . after perturbatively 
expanding both sides of Eq. (|34l) and comparing terms 
order by order, one obtains the relationship between Eu- 
lerian and Lagrangian bias parameters. To the first two 
orders [To| . 



6f(k) 
&f(ki,k 2 ) 



i + 6f(k), 

&£(ki,k a )- 
x &f(k 2 ) 4 



-F 2 (k!,k 2 )6f(k) 
[k.Li(k a )]6f(k x 



[k-ii(kr)] 

(35) 



Therefore, even if the Lagrangian bias b„ is local ini- 
tially, non-linear and non-local gravitational evolution 
J"(")(pi, • ■ • , p„) would render the Eulerian bias bf non- 
local. 



C. Continuous Galaxies Formation Model 

As briefly mentioned in the introduction, we will try to 
generalize the Eulerian bias model in two different ways. 
As the resummation technique will be discussed in next 




FIG. 6: Comparison between Lagrangian bias model (left) 
discussed in [J] and our galaxy formation model (right). The 
thin solid denotes the linear density perturbation 5q , and thick 
line represents nonlinear overdensity <5 m (n). Note that the 
right diagram has not revealed very clearly the process of 
continuous galaxy formation as shown in Eq. (1401) . 



two sections, in the following of this section, we will try to 
construct a model by combining Eulerian and Lagrangian 
pictures together. 

First, we still assume that the complicated nonlinear 
process leading to the formation of galaxies or other type 
of observable tracers can be described by a general non- 
linear functional F[S m (rj)] at the moment of galaxy cre- 
ation. In general this can be nonlocal, however, we will 
assume -F[<5 m ] to be local for simplicity. As what will 
be seen in the following, the locality assumption here is 
not as problematic as in the local Eulerian bias model 
since it only relates to the physical process during the 
moment of galaxy formation and could be easier to gen- 
eralize. Furthermore, such functional can also explicitly 
depends on the formation time, i.e. we have F[5 m (r]), 77] 
and 6„(ki... ,„,??). 

Assuming at time 77, there are Ap g (q, 77) new galax- 
ies were created during A77 at a Lagrangian-like position 
q, and then evolved to Eulerian position x where have 
been observed by a survey at time rjQ. Analogous to the 
Lagrangian bias model (Eq. |3"3"|) . this contributes to the 
Eulerian galaxies a number density p s (x, 770) 



J d 3 q Ap g (r),q) S D [x - q - * g (q, 77, n )], 



(36) 



where * fl (q, 77, 770) is the displacement field of galaxy 
from 77 to rjo- Note that Ap here only comprises galax- 
ies that are finally selected in the sample. Therefore, the 
galaxy number density at observed time 770 and comoving 
position x is a sum of all contributions that formed before 
770 and then moved to x with a displacement ^ f g (q, 77, 770), 

/no r 
dr) / d 3 q p' g (r))F[5 m (q,r)),rj\ 

x <5 D [x-q-vEr ff (q, 77,77o)]. (37) 

Here we have written the number density Ap g (q, 77) = 
( 7 ?)^[^m( ? 7): v\i where <5 m (q, 77) is the nonlinear 
evolved matter density field at 77. In principle, the dis- 
placement field of galaxies should follow the one of 



dark matter *& m , i.e. 

®gfav,yo) = by [*m(qm,?7o) - ^ m (q. m ,v)} , (38) 

where qj ra is some initial Lagrangian position of the mat- 
ter field when density perturbation is still linear. In 
the following of this paper, we will assume by = 1. In 
the simplest case where all galaxies formed at r]f, i.e. 
Pgiv) = Pg^Div — Vf), a similar form of Eq. (j3"3")l is re- 
covered. Taking the average of both sides of Eq.(f3"T|). 
one obtains the average number density p g {r]o) as a in- 
tegration of the average density over the whole galaxy 
formation history 

rvo rlo 

PgM= dri (f g (rj)(F[5 m) ij\) = dn p' g (r)) (39) 

•J TJmin •* Vmin 

Therefore, the perturbation of galaxy number density can 
then be expressed as, 

/no r 
dr) J d 3 qf( V ) e- lk t iF[ ( 5 m (q,7 7 ),7 7 ] 

x e -iK-* s (i,n,r) ) -S D (k). (40) 



where f(r/) = p g (r])/ Pgivo) denotes the average galax- 
ies formation rate, with the normalization relation 
J dr) f(rj) = 1. To first order, we have 



<^(?7o) 



dri f(v) [biD( v ) + by [D(r, ) - D(r,)}} S 



byD( Vo ) + I dr) fWDWih 



Therefore, the linear bias of galaxies evolves as 
blivo) =bq, + I dr] f(rf) [ b i - & *J 



n° Din) 



So, 
(41) 

(42) 



In the case that galaxies exactly trace cold dark matter, 
i.e. by — 1, we recover the result calculated by [ll| if 
onl y th e growing mode is considered. As first discussed 
by [la], for the model galaxies burst at a single time 77/, 
f(v) = s d(v - Vf), 



(43) 



So the gravitational growth will always debias the clus- 
tering of galaxies towards 6f = 1. 

Because of this Lagrangian- like approach adopted here, 
this framework automatically incorporates the gravita- 
tional non-locality induced after the formation of the 
galaxies. Furthermore, as will be seen in Section V., Eq. 
(|4"0")) also allows the application of the resummation tech- 
nique so that the perturbative convergence of such model 
is guaranteed. 



r (i) 

1 E ■ 





<y- 



FIG. 7: Diagrammatical demonstration of the two-point non- 



linear propagator T E (k) for Eulerian bias mode. The grey 
ellipse includes all possible ensemble averages among initial 
density perturbation 8q. 



IV. RESUMMED PERTURBATION THEORY 
OF EULERIAN BIAS MODEL 

In this section, we will show that the Eulerian bias 
model can be partially resummed with the help of ex- 

(n) 

poncntial decay of multipoint matter propagators re - 
starting from the definition of this model (Eq. |32|) . one 

can define the (n+ l)-point nonlinear propagator of 
the biased-tracer overdensity Se, 



E, f>. 



. bn (ki,--- ,k„;f7)<5_ D (k-ki... n ) 
1 / S n E E (k) 



6x(kl)- 



b„( k ") 



(44) 



Before substituting Se (Eq. I3"2")) as well as the pertur- 
bative solution of S m (Eq. [T6l IT8l) into above definition, 
one first notices that, T^ 1 ' can be expanded in terms 
of the number of nonlinear density fields S m involved, 
just as depicted in the schematic diagram of Fig.® 
for the two-point propagator. Mathematically, it means 



F 



(n, m) 
E 



, where 



F (n) _ 

1 E — 2^7, 

r B;bT- ) -fc„( kl ' ' ' ' ,k„;?7)5 D (k-ki... r 

d 3 qi... m S D (k - qi... m )6 m (qi, • • • ,q m ) 



? 1 777 I 



[S m (qi) • • • S m (q m )] ) .(45) 



For the two-point propagator (k) , let us write down 
all contributions formally. First of all, it is clear that 
the one-5 m term, i.e. the first diagram in Fig.®, equals 
the two-point propagator of the matter field itself 1^ (k) 
multiplied with the linear bias parameter bi(k). For the 
two-<5 m contribution, i.e. the second diagram, one reads 
from Eq. (gSJ) 



r£' a 2 )(k;?7) = 2 J d 3 q 12 5 D (k-q 12 
'<55 m (qi,77) 



&2(qi,q2) 



Sm(q2,f?) 



(46) 
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As discussed in Section II, the ensemble average at the proceed, we expand Eq. (|46|) explicitly 
second line can be expanded by a T— like series, which 

in) 

in this case, is the propagator for matter field T m . To 



r ^f ) ( k ; , z) = / rf3 qi2 Mk-qi 2 ) Mqi,q2) X! ("i + x ) / d 3 Pi - ni d 3 pi...„ 2 <5i)(qi - k - pi... ni ) 

ni,ii2 

x 5 D (q 2 - p!...„ 2 ) ^'S'"^ ( k ' Pi> • • • , Pm;i)-?i"l... ( i„ j (pi, • • • , Pn 2 lV) 

X (</ , ci(Pl)---0c„ 1 (Pn 1 ) <f>di(Pl) • ■ • famiPnJ) ■ 



(47) 



As we have shown, Wick's theorem ensure the de- 
composition of joint ensemble average at the last line 
into the combinations of two-point correlations. Each 
of these terms can be labeled by three indices r\,r2,t, 
where rii = 2r 1 + t, n 2 = 2r 2 + t, i.e. we classify all 
pairs into three categories: r\ pairs within the first field 
4> Ci (Pi), (1 < i < Tii); r 2 pairs within the second field 
(frdjip'j), (1 < J ' < "-2); and t pairs in between. The 
renormalization is then achieved by realizing that a pre- 
summation of r\ and r 2 gives rise to one (t + 2)-point and 
another (t + l)-point propagator, therefore 



r£a 2) (k^) =^2(U+1)\ J d 3 Pl ... n P (pi) ■ ■ ■ P (Pn) 



order n, 



& 2 (k + pi...„,-pi...„) r^+^(k, pi,-- - ,p n ;rj) 
rL n) (-Pi, - ,-Pn,v), (48) 



where we have introduced the notation 



x m V 



Pi, ' 



,Pn,v) 



An) 

- m, c\ ■ "C, 



(pi 



,p n ;ri) u c 



r£: ) (k;, ) = E^ ) {/f II 

{tij} 1 l l<i<]<n 

■■■r^(pf.. tnl ,---,P^. tnn ;v) II [ p o(p l3 )] U3 }> 



l<i<j <n 



(50) 



where the index ty denotes the number of connections 
between the z-th and the j-th density field S m (e.g. 
{tij} = {^I2,£i3,ii4, £23,^24, £34} for n = 4), so the total 
number of the indices is then n(n — l)/2. The coefficient 

<7^' ^ is then given by 



n U, 



Ki< 



n 

l<2<_7<n 



(*«)! 



(51) 



and used the definition of r„' in terms of the perturba- 
tion kernel. The coefficient n\ then comes from all possi- 
ble ways of matching n initial states with another group 
of n initial states. For this contribution, the only differ- 
ence between Eq. (|48|) and Eq. (fT8|) is the order of the mat- 
ter propagator been summed. For the RPT description of 
the matter power spectrum, both propagators share the 
same order; otherwise, it would be impossible to match 
the pair. For the same reason, the order of one propa- 
gator in Eq. (14"51) is always greater than the other by one, 
since we have to select one branch out before taking the 
average. 



Here, ti = ^2 m U m is the total number of connections 
' linked between the i-th density field S m and others, with 
(49>, y - = tji. Eq. (P) expresses the products of multino- 
mial coefficients of choosing {tn, • • • , Un} from ti for each 
density field, times all possible permutations within tij 
for pair matching. Each connection carries a momen- 
tum p 1 , 3 , (1 < I < Uj), which characterizes the i-th 
connection between i-th S rn and j-th 5 m . The integra- 
tion is taken over momenta of all possible connections 
among n different density fields. Because of the ensem- 
ble average Eq. (fT5|) , we have pj^ = — p^. 6„(k, p y ) = 
6 Jl (k + ^\p 1 \ J^i P 2 V" ,Y,iP m )> and we also used the 
shorthand notation [Po{p l ' J )] U] for Po(p\ 3 ) ■ ■ ■ Po ($>%)■ 

Given above expressions, we are then able to expand 
r^g in terms of the number of initial power spectra en- 
tering the calculation. For the tree level, it simply reads 
as 



The formula can be further generalized to arbitrary 



T^T e \k;r,)=b 1 (k) rW a (k;„;fl) 



(52) 
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FIG. 8: Nonlinear propagator T^\k) of Eulerian bias model 
up to one-loop order. The integer in front of each diagram 
indicates the number of all topologically equivalent diagrams. 



a straight path through the diagram carrying the same 




+ 



FIG. 9: Each diagram contains an infinite number of loop 
contributions. 



where we have introduced the notation 

rW(p x ... n ; 77; R) = rW ( Pl . . . p„ ; r,)W(\ Pl ... n \R), (53) 

where W(kR) is some smoothing window function. 

At one- loop level, from Eq. ((50|) . the contribution would 
be nonzero only when n < 3, since otherwise there would 
exist at least one r m with zeroth order, which would 
vanish because = (5 m ) = 0. For n = 2, £12 = 1, the 
coefficient equals 4. For n = 3, £23 = 1, ^12 = £13 = 0, the 
coefficient equals 3. So we have 



d 3 q P Q (q) 



63 (k, q, -q) 



[rW(q; m R)] 2 r« Q (k; m R) + 2 b 2 (k - q, q) 



r^ (k,-q;r ? ; J R)rW(q;r ? ;i2) 



(54) 



To make the Eq. (|50j) easier to understand, one can 
draw every contribution diagrammatically. Starting from 
the diagram representing in Fig.©, we change all 
kernels into n-point propagators. After selecting one par- 
ticular branch out, we glue the rest of the initial states 
(open circles) together. Every resultant topologically in- 
equivalent diagram represents one or several terms in 
Eq. (|50|) . Since all ordinary kernels have already been 
substituted by propagators, the ensemble average (glu- 
ing) is only performed among different density fields S m . 
In Fig.©, we show all diagrams of up to one-loop or- 
der. They correspond one-to-one to Eq. (j52} and Eq. ([5*4| . 
We also present all two-loop diagrams and equations of 
in Appendix B. It should be emphasized that, by 
substituting ordinary kernels into propagators, each dia- 
gram in fact contains an infinite number of loop contribu- 
tions at every substituting position, as shown in Fig.©. 

Numerically, the Gaussian damping of the multipoint 
matter propagator will eventually contribute to the 
improvement of the convergence in calculating nonlin- 
ear matter power spectrum. Therefore, it is important 
(n) 

to check whether T E share the similar feature. From 
diagrams representing r^(fc), we see that there exists 



momentum k at both the start and end. Along this 
path, there is one convolution bias vertex b n , as well 
as one [n + l)-point propagator rf. For n = 1 (e.g. 
the third diagram in Fig.©), this contribution recovers 

the same fc-dependence as Tm (k), rescaled by a constant 
from loop integration. When n > 1, Eq. (fl~7|) suggests a 
similar damping of k given pi ■ • • Pn-i- Meanwhile every 
propagator associates with a smoothing window function 
W(kR). Therefore one should expect that, in the large-A: 
limit, r^P will decay as a combined effect of a Gaussian 
damping of and the smoothing window function. 

Besides the tree level result of two-point propagator 
T^' tTCC \k), which equals the one of dark matter field 

times the linear bias coefficient bi(k)T$(k), the loop 
integration from higher-order calculation, e.g. the third 
diagram in Fig.©, will also change the small- A: normal- 
ization of (k) and eventually the effective large-scale 
bias of the power spectrum 6* . As we will see at the end 
of this section, T^\k) dominates the contribution of 6*, 
hence the accuracy of the 6* depends on the precision 
of estimating the r^'(fc). Due to the exponential decay 

(n) 

of matter propagators r„ , our resummed formula will 
help us in improving the convergence of b\ . 




k — v 



FIG. 10: Examples of two different contributions of three- 
point propagator, full set of diagrams is shown in the Ap- 
pendix B. 

Now, we can further derive the three-point propagator 
similarly, except that two distinct contributions have to 
be taken into account because of the second order func- 



tional derivative in the definition of T 



(2) 
E ■ 
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rg'2(k ls k 2; „) 



d 3 qi...„ S D (k - qi...„) 



&«(qi,- • • ,q« 



^ 2 ^ m (qi,?7)„ , i — i / 



+ (n - 1) 



5S m (qi,77) fc. m (q 2 ,?7)~ 



S m (q3,?y) ■ ••S m (q n ,7?) 



(55) 



As shown in Eq. ([5"5j) , the first term takes the joint aver- 
age of the products of a second derivative with n — 1 den- 
sity fields. Diagrammatically speaking, it corresponds to 
selecting two branches from a single S m , i.e. the first 
diagram in Fig. <jTUl)- In Fig. ([18]) of Appendix C, this 



category includes the first, third and sixth diagram. The 
rest of the diagrams then correspond to the second term 
of Eq. (|55p . where two branches originate from two differ- 
ent density fields. Both terms can be written in the same 
form of Eq. (|5"0)) . 



r|;)( kl ,k 2 ^) = £ / II < /:i p;' 



6n(k,p y ) 



n 



(2, n) r (ti+2)rt , li in v 

i m, ab ( k l, k 2,Pl- tll > • • ' 'Pl-dJ 



{*«} " L l<«J<ri 

x x 1 m IPl— t„i> )Pl—t nn ) ' "{tij} 1 m, a k K l > Pi— tn ' ! Pi— t lB , 



X r 



(*2+l) 



( k 2 5Pi 



21 



■«21 ' 



2n \ p(t n )/ nl 
>Pl— *2,J L m KPl—tni' 



nn \ 
• Pi— *«»J 



n 

l<i<j'<n 



(56) 



The difference between the two terms can be 
clearly seen from their orders of the matter prop- 
agators. With the same labeling system, the 
first gives rm 1+2 'rm' • • • Tm n \ while the second gives 
T^ 1+1) r^ 2+1) • ■■T { t ) - Meanwhile, the two g coefficients 
equal 



(2, n) 



-(2, n) 



-(ti+2)(ti + l) 



n 

L l<i<j<n 

n(n — 1) 



n 

Ki<n 



(tl + l)(*2 + l) 



n 



til ' ' ' tin 



n 



l<2<j<n 



(57) 



The first contribution g^ t \ gives n.(ti + 2)(ti + 1), while 
the second term gjj'.™' 1 gives n(n — l)(ti + 1)(*2 + 1)- 

At tree level, there are two diagrams: one is bi(k)T nl , 
from the single-5 m contribution; and the other is the two- 



5 m term with t±2 = for the second term of Eq. (|56p . 
rg :;° c) (P, k - p;ry; i?) = 6i (k) T%] ab (p, k - p; , /; R) 

2 



b2(k ' 9 p ' p) rw 8 (p^;fl)rg t ck-p^fl). 



(58) 



All the one-loop contributions can be found in the Ap- 
pendix C. 

Finally, the non-linear power spectrum of the field Se 
can be expressed as 

P E (k;v) = J2 n - d 3 q,i... n S D (k-q x ... n )Po(qi)---P (q n ) 



n>l 

X 



(59) 



In Fig. (ITT1) . we illustrate the numerical calculation 
of the power spectrum of Eulerian bias model at red- 
shift z = 0.5 for bias parameters b% = 1.5, &2 = 0.5 and 
f>3 = 0.2, assuming a spherical top- hat smoothing length 
R = 2Mpc//i. The solid line gives the full nonlinear 
power spectrum of Pg(fc) including both and F^ 
contributions. Each of them is then represented by the 
long dashed and the short dashed line respectively. We 
only calculate the up to one-loop order while T^' to 
the two-loop order. Compared to the dot-dashed line rep- 
resenting the linear bias contribution b\Pi,(k), one sees 
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FIG. 11: The power spectrum of Eulerian bias model at z = 
0.5 with 6i = 1.5, &2 = 0.5 and 63 = 0.2, assuming a spherical 
top-hat smoothing window function with R = 2Mpc/h. The 

solid line gives the full result Pb{U) including both and 

(2) 

T E contributions, while each of them are represented in long 
dashed as well as short dashed line respectively. Here, Ig 
is calculated up to two- loop order and T E to the one-loop 
order. The dotted line shows the linear matter power spec- 
trum Pi(fc), and the dot-dashed line illustrates the linear bias 
contribution biPi(fc). The offset between the dot-dashed and 
solid curve indicates that higher order bias parameters con- 
tribute to the large-scale amplitude of the power spectrum. 



clearly that higher order bias parameters contribute to 
the large-scale amplitude of the power spectrum. This 
can also be seen in the left panel of Fig. ([T3jl . where we 
plot the normalized two-point propagator (k, z) / 'D(z) 
for the same bias model but at various redshifts, starting 
from z = 6 to z = 0.5. At lower redshift, the variance 
of the local density field grows larger, so does the effec- 
tive linear bias b\ as shown in the figure. For the same 
reason, the exponential damping length decrease towards 
higher redshift. As expected, all of lines at large scale lie 
above the true linear bias bi, i.e. the dotted horizontal 
line. We want to remind here that the accuracy of such 
effective large scale bias calculated by our formalism has 
already been verified with the simulation in 16] for a 
logarithmic transformation. It is guaranteed by incorpo- 
rating the matter propagator, since the loop integrations, 
that contribute to b\ , are reduced and therefore the series 
expansion is regulated. 

Another distinguishing feature provided by our for- 
mula is the scale-dependent bias. In Fig. (fT21) . we give 
several examples of such scale-dependent deviation at 
baryonic acoustic oscillation sclaes compared to the lin- 
ear matter power spectrum at two different redshifts for 
various bias parameters. They are calculated by dividing 
out the no-wiggle power spectrum [l9j with appropriate 
normalizations. As can be seen in the figure, such scale 



dependence is, roughly speaking, determined by the rel- 
ative value of linear bias b\ and high-order b n . So for a 
fixed linear bias, a bigger high-order b n leads to a larger 
scale-dependent deviation, while a smaller b\ results in 
the similar trend if b n are fixed. Not surprisingly, the 
same bias model will give a larger deviation at lower red- 
shift. 



V. RESUMMED PERTURBATION THEORY OF 
CONTINUOUS GALAXY FORMATION MODEL 

Based on the construction of resummed perturbation 
theory of Eulerian bias model in the last section, we will 
further extend such model to incorporate a continuous 
galaxies formation history. As shown in Eq. (|40[) . the per- 
turbation of galaxy number density observed at 770 can 
be expressed as, 



<5 9 (k,?7o) 



S„(k, 



m ) 



no 



drj f(r))A g (\n,ri,r)o). 



(60) 



where the integrand A g (k;?y, 770) is the density contrast 
of galaxies formed at 77 and then been observed at 770 



A s (k) 




(k, 77, ?7o) - S D (k) 



I ml 



<f 5 k 1 ... n tf 3 k' 1 



n + m > 1 

x 6 D (k - ki...„ - k' 1 ...„J6„(ki, • • • ,k„;7i) 
x [<5 m (ki, 77) • • • 6 77)] [k- 1^(^,77, 770)] 
x ••• [k-* 5 (k; n , 77,77o)] (61) 

where F and e~ ik '*f are Fourier transforms of F and 
e -ik-* s respectively, and * denotes the convolution. 
Apart from the weighted average of A g with the galaxy 
formation rate /(??), this model differs from the work 
of [12| in two ways. First, the non-linear bias func- 
tional -F[<5 m ] depends on the nonlinear evolved matter 
field S m (r]) instead of initial state So- Secondly, the dis- 
placement field \& g (?7, 770) characterizes the movement of 
the galaxy from location q at time 77 to Eulerian posi- 
tion x of time 770. In the following, we will construct the 
propagator of 5 g (k, 770) instead of following the work of 
[12j | exactly. 

Starting from Eq. (|60[) . one can define the multi-point 
propagator of galaxies similarly. For the two-point prop- 
agator r^^fc), we have 



<5a g (k, 770) 
50(k') 



= rW(k, 77 )fc(k-k') 



I/O 



dllI{ll) < wy) / (62) 
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FIG. 12: Examples of scale-dependent deviation of nonlinear Eulerian bias model at scales of baryonic acoustic oscillation for 
various bias parameters at redshift z = 1 (left) and z = 3 (right). Both and are calculated up to the one- loop order. 
From the top to bottom, we adopt bias parameters: (1) b\ = 0.8, 62 = 1-5, &3 = 0.4; (2) 61 = 0.8, 62 = 1, 63 = 0.2; (3) 
61 = 1, 62 = 0.5, 63 = 0.2. 




FIG. 13: Diagrams of T g '(k) up to the two-loop order. The 
first row represents the contribution of 2 (fc) and the second 
row for r£i(ft) terms. 



Concentrating on the quantity inside the integration, and 
we define 



/ = iy(k;^o)Mk-k) (63) 



Substituting the expression of Eq. (l61[) into this defini- 
tion, one immediately obtains two different contributions 
arising from the functional derivative 

rW(k ; 77, rn) = r« (k ; v , vo) + r^ 2 (k ; n, m) 



[54>(k) 



F 



Si 

F * 



L^(k) 



tive expansion gives 



.(i) 
■g, 1 



£/dk 

n=l J 



6 n (k,ki, • • -k„) 



<S m (ki) • • • S m (k n ) 



™ — n " «^ 



[k-* ff 0O] 



(65) 



Without going into detail, one finds the resemblance be- 
tween above equation and Eq. (|46l l47|) . since both F and 
e -zk * B can wr jtten in the form of Eq. (TK)1) . There- 
fore, following the exact same procedure after Eq. (|46l) . 
the contribution I 1 ^ j can be similarly resummed using 
the T— expansion, and one simply reads from Eq. (j48|) as 



r« 1 (k; ?7 ,r7o)=^(n + l)! J d 3 Pl 



fo(Pl-n) 



x r^ +1) (k,pi,--- ,p„;?7)np(pi,--- , p»; 77, »7o) (66) 
Here the first contribution, arisen from (SF /S4>(\s)), is 

(n) 

the multipoint propagator of Eulerian biased tracer r^, 
which we have discussed in the last section, and the sec- 
ond contribution is defined as 



(»), 



ni n) ( P i, 



,p n ;v,no) 



1 I 5 n e-' ;k ^(P;^^o) 
n\ \ <50(pi) • • -<5</>(p„) 



(67) 



(64) 



Concentrating on j(k; 77, 770) first, a further perturba- 



wherc p = pi...„. When k = p, it recovers to 
the definition of multipoint matter propagators in La- 
grangian perturbation theory with galaxy displacement 
field (77,770). 
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FIG. 14: (left): Two-point propagator of Eulerian bias model T^\k) at various redshifts z/, 



0.5, 1, 2, 3, 4 and 6 from 



top to bottom. We assume the bias parameters &i = 1.5, 62 = 0.5, 63 = 0.2 and otherwise, (right): Ilj. 1 (k; z/ rm> z t s ) (Eq. 
I67p . from Zf orm = 3 and 6 evolved to z a t s = 0, 0.5, 1, 1.5, 2 and 2.5 (top to bottom). 



The physical consequence of Eq. (|66j) is the separation 
between the nonlinear gravitational evolution of galaxies 

(n) 

after their creation characterized by O^. and the more 
complicated nonlinear galaxy formation physics parame- 
terized by the unknown functional -F[<5 m ]. The final non- 
linear growth of galaxy density perturbation Tg relates 
to these two physical processes in a statistical way. Dia- 
grammatically, Starting from the second diagram in Fig. 
([6]), one can construct and Ilj: 7 ^ simply by gluing 
initial states connected to wavy line and thick solid line 
separately. We illustrate T^\(k) in the second row of 
Fig. (fT3")) up to the two-loop order. 

Similarly, the second term in Eq. 



.(i) 
9,2 



[ m 6„(ki, • • -k„;?7) 
2^ / aki-n j o m (ki) • • -d m (k„ 



E 



n=l " 
(-<)' 



1.1 



dk[... m [k-^ g (k' 1 )}-.-[k.^ g (k' m )} 



8(f>(k') 



can be resummed as 



(68) 



rW 2 (k;r ?) r ? o) = ^(n+l)! J d 3 Pl ...„ 



Po{Vl-n) 



x r^ } (k,pi,-- - ,p„;77)n^," +1) (k,pi,--- ,p n ;V,Vo), (69) 



which corresponds to the first row of Fig. ([13]). Therefore, 
up to the one-loop order, the two-point propagator of 



galaxy density perturbation T g 1 \k) equals 

If) (k; 77, 770) «r^(k; V ) +U^(k;r ] , m ) + 2 J d 3 p P (p) 



(k, p; V )U^ (p; 77, Vo ) + T { £ (p; tj)!!^ (k, p; 77, % ) 



(70) 



To proceed, we have to estimate the III™) with galaxy 
displacement field $^(77,770). To the first order, the 
Zel'dovich approximation simply gives 



(k; 77, 770) = [£>(%)- D(r,)}5 (k) 
and therefore, 

llf^fc; 77,770) » [D -£>]exp 



k 2 (D -D) 2 2 
o °* 



(71) 



,(72) 



where Do = D (770), and D — D(rf). So it is fully charac- 
terized by the gravitational growth between two different 
moments. Though higher order corrections analogous to 
the calculation in Appendix A are possible, we haven't 
explicitly calculated them here. The caveat is that the 
galaxy displacement field at q relates to ^f m at some 
initial Lagrangian position qj n since we assume 



* s (q; v, Vo) = *m(q 4 «; vo) - *m(q 2 «; v)- 



(73) 



At the lowest order q = qj n , but the coordinate trans- 
form needs to be carefully incorporated for higher order 
calculation. Furthermore, a more reliable approach is to 
directly solve the fundamental dynamical equation (Eq. 
[5]) with the potential determined by the nonlinear matter 
distribution at 77. We will leave this in the future work 
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FIG. 15: The normalized two-point propagator T g (k; Zf orm , z OBS ) / D(z b a ) of galaxies which formed at fixed formation redshift 
(left) zf orm = 3, and observed at z b s = 2, 1.5, 1 and (top to bottom); and at fixed observation redshift (right) z b 8 = 0, and 
formed at Zf orm = 0.5, 1, 2 and 3 (top to bottom). Assuming nonlinear bias model £>i = 1.5, 62 = 0.5, 63 = 0.2. Inner panel: 
the corresponding power spectra in the opposite/same (left/ right) order. 



and simply utilize the Zel'dovich approximation (Eq. l72[) 
within this paper. 

In the right panel of Fig. (fF4")l , we illustrate the Jljp (k) 
from the moment of galaxy formation Zf orm — 3 and 6 
evolve to various observational time z h s - At the large 
scale, it converges to linear result D(z ous ) — D(zf orm ), 
which relates to the linear debiasing we discussed in Eq. 
(1421) . and dampens to zero towards smaller scales with 
the speed proportional to [D(z h s ) — D(zf orm )] 2 . As can 
be seen in the figure, when galaxies have more time to 
evolve , i.e. z oos — Zf orm is large, the observed galaxy 
distribution is more influenced by gravitational evolu- 
tion. This can be easily seen from Eq. (|T0[) , since at tree 



level T^fc) 



rjg (k; Zf 



+ U ( k 1) (k;z fo 



1 Z Q bs 



When 



D(z obs ) - D(z form ) is small, Ilj^(fc; z form , z obs ) vanish 
and therefore the gravitational evolution is negligible. 

This can also be clearly seen in the left panel of Fig. 
(1151) . where we plot the normalized galaxy propagator 



Tg (k; Zf orrn , Zobs) / D(z bs) with fixed formation redshift 
Zform = 3 and observed at various z a bs- The dashed 
line shows the tree level of Eq. (fTOf and solid line corre- 

(n) 

sponds to 1-loop results. In this situation, T E is fixed, 
and as z bs decreases from 2 to 0, the large scale bias 
also decreases as expected. However, the amplitude of 
the power spectrum itself as shown in the inner panel of 
left Fig. (|15l) . is still increasing at lower redshift due to 
the linear growth D(z b s )- From the diagrams of three- 
point propagator of Eulerian bias model, at large scale 

k — » 0, (k, p) is nonzero, hence the loop integration 
also contributes to the large scale bias, although it's still 
dominated by the tree-level value. As for smaller scales, 



r^ 1 "* dampens qualitatively similar to the TL^'(k). 

In the right panel of Fig. (fl"5"j) . we illustrate the same 



quantity for fixed observation redshift z oos — 0, when 
galaxies are formed at different time with the same non- 
linear process (same -F[<5 m ] and b n ). In this case, both 
and 4 n) are changing and the evolution of large 
scale bias would in principal depends on the specific pa- 
rameters b n adopted, since the Te and life evolve with 
time oppositely. In the examples we shows, 1^ (fc) dom- 
inates the evolution, i.e. galaxies formed at lower redshift 
would be higher biased due to the stronger nonlinear ef- 
fects. 

Finally, we consider the continuous galaxy formation 
model, assuming a simple log-normal model for average 
galaxy formation rate 



ipic 



f( z ) = N g exp 



[log(l + z) - log(l + z )] 2 



2a 2 



(74) 



where N g is the normalization factor ensuring J dzf(z) = 
1. Therefore, we have two parameters characterizing the 
model: zq the peak redshift of galaxy formation, and 
a 2 the width of the redshift for galaxy formation. In 
the lower panel of Fig.lfTB"]). we give several examples of 
galaxy formation history. 

In the left-upper panel of Fig. dTBl . we present again 
the normalized propagator for various galaxy formation 
models which we shows in the lower panel of the same fig- 
ure. As discussed previously, the Tg(k; z b s ) here is sim- 
ply a time average of Tg 1 ^ (k; Zf or m, z bs) weighted with 
f(zform)- Therefore, for a narrower formation history 
(e.g. the dashed, dotted and dash-dotted lines), where 
(T 2 = 0.01, the results are similar to corresponding galaxy 
bursting models with different zq. On the other hand, 
when galaxy tends to form in a similar rate during a 
longer time (solid and long-dashed lines), e.g. a 2 = 0.1, 
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FIG. 16: The normalized two-point propagator (left) F g (k; Zf orm , z f, s ) / D(z j, s ) and corresponding power spectrum (right) 
for various models with continuous galaxy formation, (left-lower): The cumulative average galaxy formation history for models 
shown in the upper panel. 



we are more likely to observe a galaxy power spectrum 
with medium amplitude and nonlinear damping. We also 
illustrate the corresponding power spectrum on the right 
panel of Fig. (fT6|) . 



galaxies, which would introduce a effective non-locality 
in the functional -F[<5 m ]. In principle, this could be solved 
by introducing another merger contributions in the inte- 
gration, but with the cost that could further complicate 
the model. 



VI. CONCLUSION AND DISCUSSION 
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be nonlocal. For example, the Ap in Eq. ([56]) only counts 
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ambiguous due to galaxy merger. Therefore, the identi- 
fied progenitors in our formalism then include all merged 
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Appendix A: Matter Density Propagator in 
Lagrangian Perturbation Theory 



Here we consider the generalized multipoint propaga- 



tor 




When pi + • • • + p n = k, one retrieves the (n + 1)— point 
propagator TmK For n — 0, we define n(°)(k) as 



n(°)(k) 

As shown in Eq.pTl 



e- ik '*(°)) « exp 



-<l*(o)l 2 



(A2) 



nW( P ) = #)(k) 



ik * 



(P) 



(P) 



n (°)(k) TW(k,p). 



(A3) 
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We then explicitly expands the terms inside ensemble and 
average, which will be denoted as T^'Qs.) in the following 



T (1) (k, P ) 



(-i) n D r 



y 

^ n\(m-l)\ 

n=0,m=l v ' 



J d 3 q k ■ *(q) J d 3 p 



d Pl...(m-l)^(p-k- Pl . 
•••.p^.k)] [5o(p[)---5 (p' m -i)]) (A4) 



Ri 



R-2 



S (l + r 2 )(3-14r 2 + 3r 4 ) 



+ ^(r 2 - l)ln 



1 +r 



1 -r 



:(l-r 2 )(3-2r 2 +3r 4 ) 



.(i) _ 



L o, i - D(ri)kiL?\ 



(1) _ D 3 (r,)k t 
J o,3 - 2 



d 3 q Lf ) (q,-q,p)P (9) 



= D 3 (r])kikj I d'qLHq)^ (-q,p)P (9), 



(A5) 



where LPT kernel LA™)(pi, ■ • • , p n ) 



L (2) 



_p 

3P r, 2 1 

\\\ [ l ~ ^1,2] I 1 ^ Mi2, 3 ] ~ - ^1,2 



1 P 



-2/il,2M2,3M3,l] + P x T 



L< 3 ) = i[L( 3a ) +perm] 



(A6) 



Substituting the definition of L/")( P i,--- , p„) into 
Ea. (IA4l) . one can explicit carry out T"W up to the one- 
loop order, 



+ ZS^ - 1) 3 (1 



1 In 



1 - r 



(A9) 



Fig(Q} shows all the contributions up to the one-loop or- 
der. Denotes Tn]mO<-) as individual terms been summed one-loop order 
in above equation, we have 



Setting k = p in Eq. ([A7| . we get the T (1) (fc) up to 



Tiel = D(rj) 



d 3 p 

504fc 3 p ! 



Pl(p) 



6k 7 p + 5k V + 50/cV 
k ~ p 



21kp 7 + - (fc 2 - p 2 ) 3 (2fc 2 + 7p 2 ) In 



which coincide with the calculation of [2C 



k + p 
(A10) 

For three- 



point m ; (pi,p 2 ) 



ni 2) (pi ! p 2 ) 



^ 2 



(pi)50(p 2 ) 



-ik * 



■*(P) 



nW(A:) 



fc i Lf ) ( P i,p 2 ) 



2 

fc i fe^r ; (pi)ir ; (P2) 



(All) 



Assuming k = pi + P2, One recover the tree- level result 

r$(k), 



r«(k 1 ,k 2 ) = n(°)(k)F a (k 1 ,k 2 ), 



(A12) 



where i<2 is the second order Eulerian perturbation ker- 
nel. 



T, 



0.1 

(i) 
0,3 



-(i) 



D(,)^i 
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kp 



[R 1 (p) + 2R 2 (p)} 



k 2 

- -Ri(p) 



(A7) 



where we have defined a slightly different version of the 
integral functions R n (k) introduced by [4( Here we de- 
notes 



dr P Q {kr)R n . 



(A8) 



Appendix B: Two-loop Order of T^' 

In this paper, we have calculated the two-point prop- 
agator up to two-loop order. Seven non-vanishing 
contributions are depicted in Fig. tfTT]) . From these dia- 
grams, one can write down all terms explicitly 



rg; 21p) (k ;?7 )= / d 3 Pl2 p ( Pi )p ( P2 ) 



where 



(1, 21p) 



(Bl) 



/C^ 21 p) = 6 6 2 (k + pi+P2,-pi- P 2) rW(- Pl ,-p 2 ;R) 
x r< 3 ) (k, Pl ,p 2 ;fl), 
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FIG. 17: Two-loop order of the two-point nonlinear propagator. 



K. 



(1, 21p) 



= 3 6 3 (k + pi +p 2 ,-pi,-p 2 )r^) (k, Pl ,p 2 ;i?) where 



K 



(1, 21p) 



K 



(1, 21p) 



x rW(- Pi;j R)r«(-p 2 ;ii), 

4 & 3 (k-p 2 ,pi +p 2 ,-pi) T^(p 1 ,p 2 ;R) 
x r«(-p 1 ;i?)rW a (k,-p 2 ;fl) 

: 6 3 (k,pi + p 2 ,-pi-p 2 ) r$(pi,p 2 ;R) 

x r«(- Pl ,-p 2 ;i2)rW (k;i2) 



(2, Hp) 



(2, Hp) 



3 6 2 (k + q, -q)r£| a6 (k - p, p, q; i?) 
x r«(-q;i2) 



2 6 2 (k p q,p + q) rW (p,q;i2) 
xr^(k- P ,-q;R) 



(1, 21p) _ 



(1, 21p) 



6 4 (k-pi,p 2 ,-p2,-pi) rW(p 2 ;i?) 
x rW( Pl ;fl)r^ a (k,- Pl ;i?) 

6 4 (k,pi +p 2 ,-pi,-p 2 ) r^)( PljP2 ;i?) 

x rW(- Pl ;i?)rW(-p 2 ;fl)rW fl (k;i?) 



4 1,2 ' P) = g 6 5 (k,p 1 ,-p 1 ,p 2 ,-p 2 ) [r r W( Pl ;i?) 

x rW(p 2 ;i?)] 2 rW a (k ;j R) (B2) 

(2) 

Appendix C: One- loop Order of T E 

At one-loop level, five diagrams are nonzero. 



rg^Vk-p;,?) = J d 3 pP Q (q) 



.*=1 



(2, Hp) 



(CI) 



(2, Hp) 



(2, Hp) 



(2, Hp) 



2& 3 (k-p,p + q,-q) rW(q;i?) 
xfW (k-p;i?)r^ 6 (p,q;fl) 



6 3 (k,q,-q) Qb (p,k- P;J R) 



rW(q ;j R) 



7 &4(k-p,p,q,-q) r« (k-q;i?) 



rW(q;i?) "r« 6 ( P;j R) 



(C2) 



In Eq. (|C2j) . /Ci and /C 2 are two-2 m contributions from 
the first and second terms of Eq. (l55|) respectively, /C 3 
and K-4 are three- EE m contributions, and /C5 is the four- 
EE™ contribution. 
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k — p 




k — p 
— < 




k — p 




k — p 



q -q 




(2) 

FIG. 18: Three-point nonlinear propagator up to one- loop order 
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